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Abstract 

The nonlinear dynamics of collisionless reconnecting modes is investigated, in the frame- 
work of a three-dimensional gyrofluid model. This is the relevant regime of high-temperature 
plasmas, where reconnection is made possible by electron inertia and has higher growth rates 
than resistive reconnection. The presence of a strong guide field is assumed, in a background 
slab model wih Dirichlet boundary conditions in the direction of nonuniformity. Values of ion 
sound gyro-radius and electron collisionless skin depth much smaller than the current layer 
width are considered. Strong acceleration of growth is found at the onset to nonlinearity, 
while at all times the energy functional is well conserved. Nonlinear growth rates more than 
one order of magnitude higher than linear growth rates are observed when entering into the 
small-A' regime. 

1 Introduction 

Magnetic reconnection occurs in space and laboratory plasmas, when the equilibrium magnetic 
field lines break due to non ideal effects and reconnect with different magnetic topology [lj. Here, 
we consider the regime of high-temperature plasmas in the presence of a strong guide magnetic 
field, where reconnection is made possible by electron inertia and has higher growth rates than 
resistive reconnection. Such a regime is named collisionless reconnection OS]. 

Two dimensional studies of collisionless reconnection with gyrofluid models were performed in 
Ref. [1] with code REC2, in a periodic configuration with flat temperature and density equilibrium 
profiles, and more recently in Ref. [5] and Ref. [6]. Development of very narrow current layers 
was also found in the nonlinear phase of collisionless reconnection [6j [8j [9]. This suggests 
the importance of affording high spatial resolution to investigate small scale dynamics, which is 
crucial in the nonlinear phase of the magnetic island growth. Three dimensional studies with 
two-fluid model in periodic configuration were also performed in a large-A' regime [9] (where 
A' is a parameter measuring the ratio of mode wavelength and current layer width). Large- 
A' regimes, corresponding to large wavelengths, were often adopted in past investiagations for 
their apparent property of having higher growth rates, and therefore of explaining better fast 
reconnection observed in nature. Nevertheless, a comprehensive theoretical model capable of 
reproducing such observed high growth rates does not exist at present day. 

Nonlinear growth acceleration was analytically predicted for the m=l mode in tokamaks, 
corresponding to large A' regimes in slab model [10]. On the contrary, a nonlinear subexponential 
growth was predicted analytically for the collisionless tearing mode, corrisponding to small A' 
regimes in slab model Numerically however, a nonlinear growth acceleration was found in 
fluid simulations of collisionless reconnection [12] even for small A' regimes [13]. The nonlinear 
growth rate was shown to increase for smaller values of the ion sound gyro-radius p s = c s /oj c i 
- with c s = \jT e jm,i the sound speed and w c j the ion cyclotron frequency - and of the electron 
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collisionless skin depth d e = c/u pe - with c the light speed and u pe the electron plasma frequency. 
The value of the nonlinear growth rate was found up to twice the value of the linear growth rate. 

In this paper, we investigate the small-A' regime in the unexplored limit of very small p s /L±_ 
and d e /L_i_, where L± is the current width: the equilibrium characteristic scale length in the 
direction of nonuniformity. No gradients of the equilibrium pressure and temperature are set, 
therefore turbulence is not driven from the beginning of the simulations, but it is allowed to 
develop selfconsistently with the island evolution during subsequent phases of the simulations. 



2 Model 

Numerical experiments are carried out with the GEM-code (Gyrofluid ElectroMagnetic [H]). This 
code, originally written for drift wave turbulence in tokamaks, adopts an electromagnetic full- 
finite-Larmor-radius gyrofluid model of equations for electrons and ions. Temperature fluctuations 
are included in the model. The GEM-code has been benchmarked with results of Ref. [3] and 
with analytical theory in large A' regime given by Ref. [3]. The gyrofluid model is obtained by 
taking the first moments of the gyrokinetic equations [151 [T6] , with the advantage of incorporating 
finite ion gyroradius effects at arbitrary order. The zeroth and first moments of the gyrokinetic 
equation for ion and electron species (labeled here by z = i,e) are: 

cL4.ii du z » 

where the tilde symbol denotes perturbed gyrocenter densities h z and velocities u z , scalar poten- 
tial <fi and vector potential A, and temperatures T z (cf. Eqs: 34-39, 83, 86, 99-105 of Ref. |14j). 
The advective time derivative is given by d/dt = d/dt + ue • V. The time is normalized in terms 
of L±/c s , with the perpendicular scale length normalized to p s , and the parallel gradient Vn nor- 
malized in terms of L±. Therefore, the size of the contravariant magnetic unit vector component 
6 s is comparable to e -1 / 2 , where e = (qR/L±) 2 and qR and Lj_ are the scale lengths along and 
perpendicular to the equilibrium magnetic field. The parameters r z , f3 e and p z are defined by: 
t z = T z /(ZT e ), (3 e = 4:Trp e /B 2 and p z = (m z /(Zmi)). The parallel gradient is calculated in the 
direction parallel to the total (equilibrium plus perturbed) magnetic field, and the gyro-averaged 

~ 1/2 ~ ~ 

potentials are defined as <j) G = Tq' ((f)) and Q G = T ± d(/) G /dT ± . The fields are advected with the 
the E x B drift velocity = c(B x Vcj)c)/B 2 and with = c(B x V£Iq)/B 2 (where for the 
electrons we use the bare potential (ftc = 4* an d = 0). The FLR nonlinearities are represented 
both by the difference between <f>c and eft, and also by wg. The electromagnetic fields are linked 
to the sources by the polarization and induction equations: 

rfn, + #V^ + £ ^ = n e (3) 
<91ogTj_ n 

-Vily = J|| = n|| - (4) 

where the Pade approximant forms are adopted [15] : Tq = (1 — /o^V^) -1 , T^ 2 = (1 — / o 2 V 2 L /2) _1 , 
and pi = v t h,i/u c i is the ion gyro-radius (v t h,i here is the ion thermal velocity). Note that in the 
fluid limit, the gyrofluid equations are equivalent to a two-fluid model. As basic assumptions of 
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Figure 1: Time evolution of the magnetic island energy (upper figure) and growth rate (lower 
figure). The time is normalized to c s /L± and the energy is in arbitrary units. The reconnection 
growth rate is calculated as the time derivative of the logarithm of the magnetic island energy. 
Here the plasma parameters are j3 e = 5.45 • 1CP 5 , p% = p s = 0.45 d e . The nonlinear acceleration 
is found in this example in the phase t ~ 300-350, and the peak of the nonlinear growth rate is 
about double of the linear one. 



this model, we neglect fast magneto-sonic waves, and we consider only frequencies smaller than 
the ion cyclotron frequency (therefore no whistlers are present in this model). The characteristic 
scale lengths of this model are the equilibrium current width L±, the ion gyro-radius pi, the 
ion sound gyro-radius p s , and the collisionless skin depth d e . The characteristic velocities are 
the Alfven velocity v a = \J B 2 / A-KUimi and the sound speed c s = \jT e jrai^ where rii is the ion 
particle density and rrii the ion mass. 

We consider the plasma regime with A'p s = 0.75, and values of d e /p s ranging from 0.41 
to 4.1. We obtain a very thin reconnection layer with respect of the equilibrium current layer 
(d e /L± = 0.04). This requires a high spatial resolution: we use a grid with N x up to 2048 
points in the perpendicular direction and to 512 in the direction of the sheared field y. This 
also results in a low reconnecting growth rate, requiring very long runs. The parellization of the 
GEM-code into 64 processors solves the problem of the excessive time length of the runs, bringing 
the characteristic run time to a few days for the longest runs. A deuterium plasma is considered 
(p e = 0.00027) in our simulations. The value of the parallel to perpendicular scale ratio has been 
chosen as e = 18350. The ion to electron temperature ratio T{ is considered unitary: n = 1. 

The evolution of a magnetic island is studied in a Harris-pinch equilibrium, where the equi- 
librium magnetic field is a hyperbolic tangent of the nonuniformity coordinate x. We initialize 
our simulations with a periodic perturbation with wavenumber k y = 2ir/L y , where L y is the box 
size in the direction of the shear magnetic field. We add to this initial condition a white noise 
with amplitude ten times smaller than the magnetic island perturbation. The white noise has 
components in all three directions, and allows the system - whose equilibrium is two-dimensional 
- to develop structures in the third direction too. Periodic boundary conditions are imposed in 
y and in the direction of the equilibrium magnetic field s, and Dirichlet boundary conditions in 
x. The choice of the boundary conditions in x differs from those of Ref. [13], where periodic 
buondary conditions are chosen also in the x direction. The energy of the nonzonal components 
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of the fluctuating fields is measured and the growth rate is calculated as the time derivative of the 
energy (see Fig.[JJ): Gt = (dEx / dt) / (2Et) ■ Such boundary conditions allows free energy to enter 
from the borders of the simulation box by means of a nonzero gradient of the scalar potential in 
the x direction, corresponding to ExB drifts in the y directions. The consistency of our choice 
of boundary conditions has been checked by doubling the size of the box in the perpendicular 
direction x, and finding that the dynamics under investigation is not changed. 

3 Results 

The initial perturbations are let evolve in time, and different phases in the evolution of the 
magnetic island are found, defined according to the functional behavior of the growth rate. During 
an initial transient phase, fluctuations with all wave numbers interact with the magnetic island 
seed. After the transient phase, a clear linear phase is found, where the magnetic island grows 
exponentially in time and the growth rate is constant (t ~ 100 to t ~ 300 in Fig. [JJ). At the onset 
to the nonlinear phase, the magnetic island amplitude is found to grow super exponentially (t ~ 
300 to t ~ 350 in Fig. [1]). During this acceleration phase, the growth rate increases, reaches a 
peak value, whose dependence on /3 e is shown in Fig. [21 then decreases again. After the nonlinear 
acceleration phase, the island reaches a saturated amplitude and nonlinearly oscillates around 
the saturated amplitude. This corresponds to a nonlinear oscillation of the growth rate around 
the null value. During the saturation phase, the energy flows to secondary instabilities which are 
formed around the island separatrix. 

Simulations are repeated with different plasma parameters, and the value of the peak of the 
nonlinear growth rate is measured in the different cases. The nonlinear acceleration is found to 
become more prominent for higher values of Lj_/d e . In particular, we repeat the simulation as 
in the case of Fig. [TJ by varying j3 e and mantaining the same ratio of L x /p s , L y /p s and L±/p s . 
In other words, for higher j3 e , we are choosing bigger box sizes with bigger equilibrium current 
sheet's widths, whereas we keep the same value of d e . The separation in these different phases 
is found clearly in all these cases with different beta: we always find an initial transient phase, a 
linear phase, and a nonlinear acceleration followed by a saturation. The value of the peak of the 
linear and nonlinear growth rates normalized with the sound time L±/c s for these different cases 
is shown in Fig. [2j 

The linear growth rate is shown to decrease with a constant slope in the logarithmic graph for 
increasing f3 e , corresponding to a well defined power law. Increasing the value of j3 e by one order 
of magnitude, the linear growth rate decreases also by roughly one order of magnitude. This is 
consistent with analytical predictions, describing the reconnection process as occurring in a thin 
"inertial" layer with relative width comparable with d e /L±. The smaller the value of d e /L±, the 
less effective is reconnection. The peak value of the nonlinear growth rate is a measure of how 
fast the reconnection process occurs during the burst phase. In Fig. [2] we show that, contrary to 
the linear growth rate, the nonlinear growth rate dependence on j3 e becomes weaker and weaker 
for higher f3 e . This means that the nonlinear acceleration becomes more relevant the higher is the 
value of beta. This also means that reconnection process has a tendency to become independent 
on the microscopic scale size. Increasing the value of f3 e by two orders of magnitude, the non- 
linear growth rate decreases by only one order of magnitude. For computational constraints, the 
maximum value of (3 e that we can investigate at present time is j3 e = 1.6 • 10~ 3 . For this value of 
/3 e , the nonlinear growth rate is 'JnlL±/c s = 0.009, to be compared with the linear growth rate 
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Figure 2: Nonlinear growth rate (upper points, blue) and linear growth rate (lower points, black) 
normalized with the sound time L±/c s , versus the electron beta f3 e . In these different cases, d e 
is the same and the other characteristic lengths are rescaled keeping constant their ratio with p s 
(therefore the sound time does not depend on f} e ). Explosive reconnection is observed for values 
of beta of the order of j3 e ~ 10~ 3 (corresponding to small- A' regime), where the nonlinear growth 
rate is more than one order of magnitude higher than the linear growth rate. 



jlL±/c s = 0.0008. We call such a regime explosive reconnection. 

The evolution of the fields during the nonlinear acceleration phase has been investigated, and 
some markers have been recognized as characteristic of this particular phase. The analysis of their 
time evolution shows that they evolve with the same characteristic time scales of the magnetic 
island, with no higher frequency component. At the O-point, secondary vorticity structures are 
found by taking a cut along the line y = 0, s = 0. The vorticity structure during the whole linear 
evolution has definited sign on each side of the plane x = 0, consistent with a uniform inflating of 
the magnetic island. On the other hand, we find an inversion of sign inside the magnetic island 
starting with the end of the linear phase. In the case of /3 e = 5.45 • 10 -5 , pi = p s = 0A5d e , the 
linear vorticity structure is doubled and reproduced inside the separatrix during the nonlinear 
acceleration phase, while the growth rate passes from the linear value to double of the linear 
value. At the X-point position, the formation of a Sweet-Parker like current sheet is found, 
corresponging to a structure with two Y-points. A successive recovering of the X-point structure 
is found at the beginning of the saturation phase. We find that current layers form during the 
linear phase along the whole separatrix, and their amplitude and steepness increase strongly 
during the nonlinear acceleration phase. The end of the acceleration phase is characterized by a 
saturation of the current spike amplitude, and by the formation of new current structures inside 
the magnetic island. Acceleration in the magnetic island growth are also found to be linked to the 
presence of turbulence around the X-point, having the possible effect of multiplying the number 
of X-points and consequently the rate of reconnection (similarly to results shown in Ref. |17|). 

The ion finite Larmor radius (FLR) effects are studied by varying the ion to electron tempera- 
ture ratio Tj. We find that in the linear phase the FLR effects yields a slightly faster growth rate, 
consistently with the analytical predictions [3]. The same is valid also in the nonlinear phase. 
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Nevertheless, no qualitative difference is found in the case of = and r» = 1. The reason is 
that ion moments don't enter the physics for J'-driven processes at low-beta regimes [4j, namely 
when k\\VA S> fe||C s , and for small A'p s . On the other hand, FLR effects were found to be more 
relevant for 2D models in large A' regimes [3]. 

4 Energy consistency 

The check of energy consistency is of crucial importance in reconnection simulations. In fact, 
reconnection processes can occur in numerical experiments also because of artificial - nonphysical 
- resistivity, such as numerical dissipation. In general, an artificial resistivity causes that the 
terms of Ohm's law calculated separately don't balance exactly in the subsequent time steps. 
Therefore, a numerical dissipation can act as a virtual resistivity and the result is a fictitious 
reconnection, whose growth rate can be in some cases higher than the physical one, depending on 
the chosen plasma parameters. In nonlinear simulations, smaller and smaller scales are created 
due to nonlinear coupling of macroscopic modes, and it's important to avoid the accumulation 
of energy at the scales of the numerical grid. To this aim, hyperviscosity is added to our model 
equations, to filter small scales physics out, just above the grid size. 

The contribution to the reconnection growth rate of subgrid dissipation, given by hypervis- 
cosity and numerical dissipation, is calculated with a specific diagnostic in the GEM-code. This 
error growth rate Ge is defined by: 

Ge = Gt — Gso + G,si (5) 

Here Gt is the reconnection growth rate, calculated as the time derivative of the total energy 
of the system: Gt = {<IEt / dt) / (2Et) ■ The total energy is the volume integral of the sum of 
the energy densities Ue, Ut, U v and U m , being respectively the ExB energy, the thermal state 
variable energy, the flux variable energy and the magnetic energy [TJ]. Gs is the sum of the 
source growth rates: Gj = /i e up£:V Jo, G a = —4>V\\J and G r = p e V\\J. Gsi is the sum of the 
absolute value of the damping rates, namely Landau damping of ion and electrons (where the 
electron Landau damping is the relevant one for our regime). The contribute of the numerical 
error Ge is checked to be always negative during the whole run. In this case, the contribute of 
the subgrid dissipation to the reconnection rate is null. 

5 Conclusions 

In summary, we have presented the results of simulations on collisionless reconnection in a 3D 
Harris-pinch equilibrium, with flat density and temperature initial profiles. The gyrofluid code 
GEM has been used, and single helicity reconnecting modes have been considered in a small 
A' plasma regime. The limit of small p s /L±, d e /L± has been investigated. The acceleration 
phase during the nonlinear growth has been analyzed and the nonlinear growth rate peak scaling 
has been provided, with respect to the equilibrium plasma parameters. Numerical evidence 
of explosive reconnection has been found for j3 e ~ 10 -3 . The underlying mechanism of our 
observation is thought to be linked with the presence of the high velocity fields. In fact, a 
collisionless regime allows the formation of strong velocity fields during the nonlinear acceleration 
phase. These velocity fields are retained to be responsible for the splitting of the X-point into two 
Y-points, and as a consequence, the explosive regime is entered. The presence of an acceleration 
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in the growth of the magnetic island, could be important in explaining the fast growth of magnetic 
islands as seen in nature. 

Simulations have been repeated with double or half amplitude of the initial magnetic island 
seed, and no difference in the value of the nonlinear growth rate has been found. Some markers of 
the acceleration phase have been pointed out. At the O-point, the formation of smaller vorticity 
structures inside the magnetic island has been found. At the X-point position, the X-point has 
been shown to split into two Y-points during the acceleration phase, and then to become again 
an X-point during the saturation phase. Along the whole separatrix, current spikes at scales 
narrower than the collisionless skin depths have been found. The importance of these effects 
in causing the reconnection acceleration is under investigation, and a detailed analysis will be 
published in a dedicated paper [18]. 

It is not quite understood whether the three-dimensionality of the configuration is a key ingre- 
dient for nonlinear acceleration. In fact, 2D reconnection is known to be governed by conservation 
laws: being cj) and A functions of the perpendicular variables x and y only, then the current J 
and the vorticity O are directed in the direction parallel to the guide field, and there is no inter- 
action between neighbour helicities. On the other hand, in 3D simulations neighbour helicities 
can interact - for instance J can interact with Vp and V(/> - and the nonlinear phase is predicted 
to have a qualitatively different evolution. 
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